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Abstract. We present a domain decomposition method (DDM) devoted to the iter- 
ative solution of time-harmonic electromagnetic scattering problems, involving large 
and resonant cavities. This DDM uses the electric field integral equation (EFIE) for 
the solution of Maxwell problems in both interior and exterior subdomains, and we 
propose a simple preconditioner for the global method, based on the single layer op- 
erator restricted to the fictitious interface between the two subdomains. 
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1. Introduction 

Solving scattering Maxwell problems in harmonic regime can be achieved with various 
methods, among which integral equations (which lead to the so-called boundary element 
methods) have proven their efficiency. Their main advantage is that they allow to 
replace a problem posed on the whole space by an equation posed on the surface of 
the scattering obstacle, reducing a three-dimensional problem to a bi- dimensional one. 
With the development of such methods, several difficulties arose successively : 

• These formulations lead classically to linear systems involving dense matrices 
(in contrast with finite element methods, for instance). Several methods among 
which the most famous is probably the FMM (Fast Multipole Method) [28] , |29j 
have been used to circumvent this difficulty. 

• There might exist irregular frequencies for which the problem is ill posed [26J. 
This is typically the case for the so-called EFIE and MFIE formulations. Other 
types of formulations (e.g. the CFIE) are instead well-posed at any frequency 

• The desire to deal with high frequency problems imposes to use fine discretiza- 
tions of the equations and consequently to solve large linear systems. This 
prevents the use of direct solvers, and one usually employs iterative methods. 
On the one hand, this needs a fast matrix-vector multiplication (which is often 
realized through the FMM), while on the other hand iterative methods become 
sensible to the condition number of the system. It has been shown that the un- 
derlying systems arising from integral equations are usually badly conditioned 
and there is a need to develop preconditioning strategies in order to accelerate 
the convergence of the iterative solver [TU], [52], [30], [3T]. For instance, the 
so-called GCSIE methodology has been developed which turns out to be partic- 
ularly efficient in the case where the object has no cavities and no singularities, 
by building intrinsically well conditioned integral equations p], [6J, [H], [23] , 

mi- 
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Nevertheless, when facing realistic problems, one has to treat large objects with com- 
plex geometries and new problems are encountered. In this paper, we address the impor- 
tant issue of resonant cavities, motivating the use of a domain decomposition method. 
Indeed, this is a particularly crucial problem in stealth applications as one needs to take 
into account the existence of large and resonant cavities, such as air intakes, or cockpits 
for aircrafts. In classical numerical computations of radar cross sections, these cavities 
are usually closed in order to avoid the poor convergence of the algorithms jl] giving 
unrealistic results. 

In this paper, we explore a new strategy to deal with this problem. Indeed, we 
intend to use a domain decomposition method (DDM) in order to split the exterior 
domain into two subdomains one of which being the cavity. The aim is to decouple the 
exterior problem (without any cavity) from the problem with boundaries (the cavity 
itself). This introduces an artificial interface X between these subdomains and a new 
coupling problem posed on E (Fig. [TJ. For simplicity, we here use on each subdomain 
the EFIE to solve the corresponding subproblems and to couple the solutions on S. 
This naive DDM algorithm turns out to converge badly. In a latter part we propose a 
preconditioning technique to accelerate significantly the solution of the DDM. 

Historically, the first domain decomposition methods for Helmholtz or Maxwell prob- 
lems were applied using a finite element method (FEM) in the interior bounded sub- 
domains and a boundary element method (BEM) in the exterior unbounded domain. 
For instance, Hiptmair considers FEM-BEM methods, first applied to acoustic problems 
|20| and then to electromagnetic problems |21j . For Helmholtz transmission problems, 
domain decomposition methods have been used by Balin, Bendali and Collino [1] to 
specifically treat the case of an electrically deep cavity, and an integral preconditioner 
using the Calderon formulas has been developed by Antoine and Boubendir [3j. For 
Maxwell transmission problems, Balin, Bendali and Millot |5j on the one hand, Collino 
and Millot on the other hand |12j . propose algebraic preconditioners which use 
overlapping or nonoverlapping domain decomposition techniques. In iterative domain 
decomposition techniques, which are split into overlapping and nonoverlapping DDM, 
the subdomains classically exchange Dirichlet or Neumann data. A substantial improve- 
ment using absorbing boundary conditions is made by Despres [15], |16j . The Schwarz 
method, originally used with Dirichlet or Neumann conditions for overlapping domains, 
is then adapted by Gander, Halpern and Nataf [19], to nonoverlapping subdomains with 
more general conditions, of Robin type. The resulting algorithm converges with a high 
convergence rate for the wave equation in dimension 1. Gander, Halpern and Magoules 
[18] optimize the method by taking more general conditions, for the Helmholtz problem 
in dimension 2. Eventually, Dolean, Gander and Gerardo-Giorda [T7] adapt it to obtain 
a Schwarz optimized method for the harmonic Maxwell problem in dimension 3. 
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We present here a nonoverlapping domain decomposition method. This DDM couples 
the subdomains through the help of an operator, instead of transmitting at each itera- 
tion the appropriate conditions from one subdomain to another. We use only integral 
equations to solve the boundary value problems in the subdomains. In particular, the 
interior problem is treated with the help of an integral equation, instead of a more clas- 
sical finite element method. We are not aware of the use of such techniques for solving 
Maxwell equations with the DDM in the context of integral equations in the literature. 



Figure 1. The scattering problem (left) and the decomposition of 
domain f2 (right). 

The paper is organized as follows. The scattering problem is first described and in a 
second part, we present the domain decomposition method and the condensed problem 
on the interface. The third part gives a quick overview on classical integral equations 
methods, and especially of the one we use here, namely the EFIE (Electric Field Integral 
Equation). The Dirichlet-to- Neumann map of the interface £ plays a very important 
role that we describe carefully in the fourth part and this enables us to present a simple 
analytic preconditioner for the employed DDM, in the fifth part. A validation of the 
method is presented using pseudo-differential calculus. Eventually, the sixth part gives 
some numerical results. Substantial improvements are shown validating the approach. 

2. The boundary value problem : assumptions and notation 

We consider a compact set D with a smooth boundary Tp. We are particularly in- 
terested in the case where the set D contains a large cavity, as illustrated on Fig. [T] 
We assume that the open exterior domain f2 = M 3 \ D is connected. 

Our purpose is to solve the harmonic Maxwell problem when D stands for a scattering 
metallic object |13j . |26j . Waves propagate with constant wave number k in the exterior 
unbounded domain £1. The electric field E is a vector- valued function which satisfies 
the harmonic Maxwell equation 



while the related magnetic field is given by 

(2) H = ^VxE. 

ik 

An electric field is said to be radiating if it satisfies the well-known Sommerfeld radiation 
condition 




(1) 



Vx (VxE) - k 2 E = inO 



(3) 
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Let E; nc be an incident electric field, the electric field E scattered by the obstacle D 
is the electric radiating field satisfying the boundary condition 7dE = — 7£>Ei nc on Fd, 
where 7/5 = npx is the metallic trace onTu, rip being the unit normal outward to D. 
In other words, the field E is solution of the following problem 

{Vx(VxE)-fc 2 E = in Q, 

7i)E = -7i)E inc on T D , 

lim (x x (VxE) + z/c||x||e) = 0, 
||x||— >oo V / 

usually named as the perfect electric conductor (PEC) problem. 

3. Notation for the domain decomposition method 

Domain decomposition methods rely on splitting the computational domain into sev- 
eral subdomains. We present hereafter the application of the method for our case when 
f2 is decomposed into two subdomains f2 + and Namely, we introduce an artificial 
boundary surface S, which splits £1 into an interior bounded domain Q~ and an exterior 
unbounded domain Q + (FlG. flj). We denote by = n dfl^, in such a way that 
the boundary of Cl^ is U S. We call the inward unit normal to £1^. The notation 
(Jq = n^x and cr^ = j^n^ x (Vx) stand for the classical electric and magnetic traces 
on E. 

We introduce the short-cut field E sc , which is the radiating electric field defined on 
Q + , having a tangential trace on U S, and such that 7dE sc = — 7£>Ei nc on and 
fj^Esc = — (T^Ei nc on r (Fig. [2]). In other words, E sc is the field scattered by the object 
when interface £ becomes metallic. 




Figure 2. Short-cut field E sc . 

We denote by E^ c the restriction of the incident field Ej nc to the domain f^. We 
look for the scattered field solution of the PEC problem Q under the form E~ — Ej^ c 
inside and E + + E sc inside f2 + , where E + and E~ respectively belong to spaces of 
admissible waves W + and W~ . 

More precisely, the space W~ is the set of all electric fields E~ which are defined on 
Q~ , have a tangential trace on U S, and satisfy 7dE~ = on r^>. Correspondingly, 
W + is the space of all radiating electric fields E + which are defined on Q + , have a 
tangential trace on U S, and satisfy 7_dE + = on Tj. Since the subdomain 0~ is 
bounded, the radiation condition is not required for the fields in W~ . 



A SIMPLE PRECONDITIONED DDM FOR ELECTROMAGNETIC SCATTERING PROBLEMS 5 



The total electric fields therefore have the expression 

E ,~. ■ E~ in Q- 



E + E inc + E S c infi+ 

whereas the total magnetic fields (computed from the electric fields with Q) have the 
expression 



H 



tot 



H inO-, 
H+ + H+ C + H sc infi+. 




= E + +E+ C + E SC 
= H++H+ C + H SC 

Figure 3. Total fields E tot and H tot . 

Since E is an artificial boundary, the total fields E tot and H tot are continuous across 
E, and the problem Q becomes the transmission problem 
(5) 

fn+ x E- = n+ x ( E+ + E+ r + E SC Y 
Find (E + ,E~) G (W + ,W~), < V J \ on E. 

n+ x H = n+ x H+ + H+ + H s ^ 



Notice that by construction the short-cut field E sc verifies n + x Ej£ c + n + x E sc = 
on E, and thus, defining the right hand side current 



u rhs = n + x EEL + n+ x H sc on E, 



equation ^ rewrites as 

(6) j E tan = E t + an> Qn 

I -n x H = n+ x H+ + u rhs , 

The preceding system, which expresses the DDM, will be solved using integral equa- 
tions inside each subdomain. We recall these integral equations methods in the next 
section. 

4. Integral equations 

Integral equations methods are commonly used to solve electromagnetic scattering 
problems. We hereafter give a short overview of the construction of these methods. We 
first recall the definitions of the single and double layer potentials, as well as the funda- 
mental Stratton-Chu formula, before describing the principle of those integral equations. 

Let Dq be a compact and connected subset of M 3 with a smooth boundary Fq , defin- 
ing two open and connected domains, the interior bounded domain J7q and the exterior 
unbounded domain Qq. We denote by n the unit outward normal to Fq and by 
the tangential trace on Fq from domain Qq . The PEC problem can be formulated as 
follows: find the electric radiating field E defined on Qq and satisfying the boundary 
condition 7E = uo, where 7 = n x 7^ is a trace on Fq, and uo = — n x Ej nc is a given 
current depending on an incident field. 
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The classical vector potential Q maps a tangential vector-field u E 2?^ (To) to the 
vector-field defined on Qq and S7q by 

1 f e ik\\x-y\\ 

(7) Gu(x) = - — / . My) dy, 

47r Jr„ If - y\\ 

where ||.|| denotes the euclidean norm on IR 3 . Then we define the single layer potential 
T and the double layer potential K. by 

(8) T=|vx(Vx5) and £ = Vx£. 

ik 

The electromagnetic potentials satisfy the following important property: given a current 
u on To, the fields Tu and )Cu are automatically solutions of Maxwell equation and 
the radiation condition ([3]), [26]. The boundary operators T and K are obtained from 
the electromagnetic potentials and are defined by 

(9) n x T = n x 7t(T) and n x K = n x jt(K,) + Id /2. 

It turns out that T and K are pseudo-differential operators respectively of order +1 and 

-i ®, m, m- 

The Stratton-Chu formulas [T3], use the single and double layer potentials to 
express an electric radiating field E and the related magnetic field H in terms of their 
boundary traces. 

(10) E = T(n x H) - JC(n x E) and H = -K(n x H) - T(n x E). 

These formulas, also known as representation theorem, are the foundation of integral 
equations, as we shall see now. 

The incident field does not satisfy the radiation condition and therefore the represen- 
tation theorem (10) does not apply to E; nc and Hi nc . Instead, because E; nc and H; nc 
are continuous on the whole space M 3 , their traces have no jump across To and one can 
show that 

(11) = T(n x H inc ) - £(n x E inc ) and = -/C(n x H inc ) - T(n x E inc ). 



Summing up (jlOl) and (11), combined with the PEC boundary condition, we obtain the 



mc ■ 



EFIE and MFIE equations, 

(12) EFIE : Tu = -Eg£, MFIE : |nxK+ildJu = nxH 

where E|^° is the tangential component of Ei nc and the unknown u is equal to n x 
7T (H + H inc ). 

Unfortunately, the EFIE and the MFIE are well-known to be ill posed at resonant 
frequencies [26]. Their linear combination, weighted by an arbitrary parameter a G]0, 1[, 
yields the CFIE, which instead is well posed at any frequency [23] , [7], 

(13) CFIE : (1 - a)Tu + Q^nxK + ^Idju = -(l- a)EgJ + an x H inc . 

In what follows and for the sake of simplicity, since we are mainly interested in the 
interface problem on S, we only concentrate on the EFIE for solving the electromagnetic 
problems inside the subdomains. 
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5. Admittance operators and the DDM 

The structure of the problem ^ naturally leads to introduce the so-called admittance 
operator^ on X 

A± : E&, -> n± x H± 

where H ± = j^VxE ± , and E ± £ W ± solves 7tE ± = E^ n . Notice that although the 
input and output data of A^ are defined only on X, the admittance operators are highly 
non-local and depend on the whole geometry of the domains £1^. 



We remark that 



A+ = R ± A ± P ± , 



where : 2?'(X) — > D'iT^ U X) extends by 0, on dQ^, data defined on X, while con- 
versely, : P'(r^ U X) — > 2?'(X) restricts to X data defined on dtl^. Here are 
the admittance operators of dfl^ which map currents E^ n defined on the whole closed 
boundaries dfir^ to their magnetic traces n 1 * 1 x H 1 * 1 . 

The system (|6|) can then be expressed in terms of the admittance operators A^ 



(14) 



E~ = E + 

-■-'tan -■-'tan ' 

— A s E tan = A^E^jj + u rns , 



which eventually reduces to 

(15) (A++A E )E tan = -u rhs , 

with E tan = E+ n = E~ an . 



Equation (15) is at the heart of our domain decomposition method. We explain below 
how the admittance operators A^ and A^ are numerically computed, while Section [6] is 
devoted to the preconditioning of the subsequent linear system. 

The admittance operators A^ can be naturally obtained by solving an integral equa- 
tion, which involves four electromagnetic potentials described below. The main differ- 
ence between our particular case and the classical theory is that the domains £l + and 
0~ are not complementary one to another. The boundaries of these domains are there- 
fore distinct, although they share the same interface X. Consequently, the convolution 
operators with the Green kernel related to the exterior and the interior electromagnetic 
potentials, that we next introduce, are not defined on the same surfaces. 



Similarly to the potentiel Q defined by ([7]), we define the vector potentials which 
map tangential vector-fields £ T>' T (T^ U X) to the vector-fields defined on by 



p iK\\x — y\\ 

47r 7r±us \\x - y\\ 

As before, the potentials are used to define the single layer potentials and the 
double layer potentials IC^ as follows, 

= --Vx (Vx^) and 1C± = VxQ ± , 
ik 



^Such operators are also classically called Dirichlet-to-Neumann or Steklov-Poincare operators. 
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while 

x T± = x ^(T*) and n ± X = X 7r(/C ± ) + Id /2, 

where 7^7 stand for the tangential traces on d£l + and d£l~ . As previously, as long as the 
boundaries U £ are smooth, pseudo-differential operators and are of order 
+1 and —1 respectively. 



Given an electric field E 1 * 1 E l/F 1 * 1 and its magnetic counterpart H 1 * 1 , we define the 

±T7± _ „± ^ TT>± „^ ^±T7>± _ „± ^ I I ™ .£)r>± _ T± 



following electromagnetic traces on the boundary U £ 



(16) a^E ± = n* x E* and erf E = n* x on dtt* = U S. 



Using the representation theorem (10) in the domain J2 , we obtain the expression of 
any electric field E 1 * 1 in in terms of its boundary traces on U S. 



(17) VE± G W , E± = 7* (of E±) - /C ± (a±E ± ). 

Now, for a current E P^(E) defined on the fictitious interface S, we have A^Uq = 
cr^E 1 * 1 ) where E 1 * 1 E is such that E^ n = P^Uq . Several integral formulations 
can be used to compute effectively A^Uq . As an example, we shall see hereafter that 
A|u^ = R^-u* where T^u*) = (Jld+K^x) (P ± u^ : ). 

Indeed, restricting ourselves to the exterior subdomain, our first goal is to find E + E 
W + such that E^ = P + u^j~. Applying the trace o-q to the Stratton-Chu formula (17) 
leads to 

(n+ x T+)(cr+ E+) = Q Id +n+ x K + ^j (o+E+) on dtt + . 
Taking the cross product of the previous equation with — n + yields 

T+(n+ x H+) = Qld+K+n+x^ (E+J on dtt + . 

This problem has the form of an electric field integral equation (EFIE): 
(18) 

Find the current u + , such that T + (u + ) = Qld+K+n+x^j (P + u^) on dn + . 

The restriction to £ of the solution u + of (18) is the trace n + x H + we look for, and 
we therefore have 



A+u+ = P + tf 



as claimed. 



The cavity (O ) is treated similarly, with the restriction that A s is well defined. This 
is the case when k 2 is not an eigenvalue for the interior Maxwell problem. 

Remark 5.1. Since this method is based on an EFIE formulation, it applies to the 
computation of the admittance A^ in both subdomains, but with the following caveat: 
metallic problems having irregular frequencies fll3]. |26j ). the EFIE is ill-posed at fre- 
quencies close to these resonances. Despite this drawback, the EFIE is still widely used 
for it is one of the methods which give the most accurate results. 
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6. PRECONDITIONING THE DDM 

As we shall see in Section[7J the equation (15) unfortunately leads after discretization 
to an ill-conditioned linear system. We therefore propose a simple preconditioner in 
order to obtain a tractable DDM which improves the convergence rate. 

The purpose of this section is to introduce a theoretical framework which suggests 
that the preconditioned equation is well posed. Unfortunately, these theoretical results 
only apply so far in an ideal setting which is not satisfied in practical situations. They 
should therefore be regarded as a heuristic and hopefully as the foundation of future 
more general results. 

6.1. A preliminary lemma. Let -Do be a compact subset of M 3 with a smooth bound- 
ary To, defining two open domains : the interior domain (the interior of Dq) and the 
exterior domain Uq = M 3 \ Dq. We define as before the admittance operators related 
to the boundary To, Aq for the interior domain and Aq for the exterior domain. Also, 
the single layer potential T and its tangential trace T are respectively given by Q and 
([9]) for the boundary IV 

The idea for preconditioning our method is based on the following lemma. 

Lemma 6.1. If k 2 is not an eigenvalue for the interior Maxwell problem, then Aq is 
well defined and we have 

(19) (A++A )T = Id, 

(20) T(A+ + A ) = Id. 

Proof. Let u 6 D' t (Tq) be a current on Tq. We define E = 7~u on J7 + and Q~ . Then by 
continuity of the potential T across Tq, 

E tL = E tan = Tu - 

Thus 

n+xH+ + n-xH- = A+E+n + A E t ~ n = (A+ + A )(Tu) 

= (cr+T - a^T)u = u, 
since the Neumann gap of the single layer potential is the identity |26|. We recall that 



af 1 is the electromagnetic trace defined in (16). We have proven (19). 



For (20), since k is not an eigenvalue for the interior Maxwell problem, the operator 
T is bijective and there exists v E T>' t {Tq) such that u = Tv. Defining E = 7~v on fi + 
and 0~ leads to 



E t + an = E t " an = Tv = u. 



Consequently, 

T(A+ + A )u = T(cj+ - <rf )Tv = Tv = u, 
since the Neumann gap of the single layer potential is the identity. □ 



This lemma suggests to precondition the equation ( 15 ) by the operator Ts defined 
by the operator T restricted to the interface X. 



Definition 6.1 (Operator Ts). 

We denote by F the explicit kernel of the single layer potential T , which can be computed 
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from and Q). The operator Ts : T>' T (Y1) — > T> ; T (E) is defined as the convolution 
operator, restricted to E 

T s u(x) = J F(x - y)u(y) dy. 
Notice that Ts does not depend on T^. 

6.2. A preconditioner for the DDM. We aim at proving that the operators T^(A^ + 
A^) and (A^ + A^)Ts are compact perturbations of the identity, using arguments of 
pseudo-differential theory. This unfortunately restricts our results to smooth boundaries 
dtl^ which in turn implies that <9f2 is not smooth in general (see Fig. [I]). 

Therefore, we consider a simplified setting in which we assume that the boundary 
is not smooth but such that both boundaries T J U £ and U S are of C°° regularity. 
For instance, in dimension 2, this implies the existence of two cusps (FlG.El on the left ). 




Figure 4. The setting for the theorem (left) and the real case (right). 



Let So be a compact subset of £, and let xo : £ — > [0, 1] be a C°° cut-off function, 
supported in the interior of E, and such that xo = 1 on So- We denote by Ts = 
^>±rp±p± anc j the operator 

(21) fs = XoTs- 

We prove that Ts is a good left preconditioner when applied to functions supported 
on Eq. 
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Theorem 6.1 (A preconditioner for the DDM). 

Let Hf, )Q (S) = {u G H^(E), such that u = on E \ E }. For u G Hf, (S), we 
/lave 

(22) fs(A+ + A s )u = u + v, 

where v G H^ 1 (S). More precisely, Ts(A^ + A^) is a compact perturbation of the 
identity in H^ (E). 



(n+ x T+)(-A+n + x) - ( n+ x K+ - - Id ) = Id . 



Proof. The representation theorem (17) applied to the smooth boundary U E yields 

1 
2 

Therefore, 

(23) T+A+ = i Id +K + n+ x . 

Since the boundary U E is smooth, the operator K + n + x is of order —1. Our goal is 
to extend this result to the operator T^A^. We have 

fsA+ = ( Xo R + T + P + ) (F+A+P+) 

= xoF + (T + A+)P + + xoF + T+ (P + R + A + - A+) P+. 

It is clear to see that for u G Hf. (£), v := (P+R+A+ - A+) P+u = (P+R+ - 
Id)A + P + u vanishes on E. Since T + is a convolution operator with a kernel F(x,y) 
which is C°° for x ^ y, we have xoF + T + v G HS?(E). This shows that the operator 

D+ = X0-R+T+ (P+F+A+ - A+) P + 

is of order — oo. Note that we have used the fact that the support of xo is included in 
the interior of E. 



On the other hand, (|23|) leads to 



XoF + (T + A + )P+ = X oR + ( ^Id+K+n+x J P + = X o^ Id+x F + (K + n + x)P+. 

Notice that xo - Id = -IdinHy (E) and that xoF + (K + n + x)P + is a pseudo-differential 
operator of order —1 in H^ (E), and is therefore compact. Remark that we have used 
the fact that P + (Hf, (E)) C H s T (dQ+). 

Having the same results for the interior case, we obtain on H^ (E) 

Ts(Af) = hd+ Xo R ± (K ± n ± x)P ± + D ± , 

where xoF^K^n 111 x)P =t and are pseudo-differential operators of order —1 and — oo 
respectively. This concludes the proof. □ 

In the real case, is smooth and thus the boundaries U E are both lipschitzian 
but not of C 1 regularity. Therefore the operator is no longer a compact operator in 
U S T (r| U E) . To study this first possibility is to come back to the case of the 

theorem by distorting the boundaries T^UE in new boundaries T^UE such that these are 
C°° (see Fig. [4j on the right), and by introducing the cut-off function xo- The theoretical 
analysis of these two approximations (the change of boundaries and the multiplication 
by the smooth function xo) is n °t straightforward. A more direct approach would be 
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to extend the theory of integral equations on surfaces with singularities. Such a theory 
was developed for bi-dimensional Helmholtz problems in |25j , but its extension to three- 
dimensional Maxwell problems remains to be done. 



7. Numerical results 

In this section, we first describe the numerical discretization chosen for the admit- 
tance operators and A^ . We then explain the discretization of the preconditioning 



by the operator of the equation (15), which couples the subdomains in our domain 
decomposition method. 

We want to solve a numerical discretization of 

(A^+A^)u = uo on S. 

We recall that A|v = i? ± v ± where v ± is solution of the EFIE: T ± (v ± ) = (± Id+K^i^x) (P ± v ) 
To describe the action of the operator A^, one needs to discretize the EFIE. 

We denote by a family of triangulations of 80,^ = T^US such that := T^nTr 
is a family of triangulations of E. The space of ffdiv-conforming Rao-Wilton-Glisson fi- 
nite elements on is denoted by Xh, and we denote by (</?i)i<i<Ar its basis functions. 
Similarly, the notation stand for the spaces of RWG finite elements on T^, asso- 
ciated with basis functions (tpf 1 )i<j<jv±) where > N and where we assume that 
Mi < N, i\)f = (fi. These assumptions will allow us below to apply the preconditioner 
to a vector in the space X^. 

Let us describe the numerical computation of operator A^ applied to a vector vo = 
Y^i=i v o,i l Pi of Xfr. We first extend vo to a vector of X£ defined by 

(N N+ \ 

i=l i=N+l J 

and we denote by 

Vq = Oo,i, • • -,vo,n,0, • • • ,0) T 
the vector of its components in the basis of X^ . 

Given an operator A and a space X of RWG functions associated with a triangulation 
To, we denote by [A]x its Galerkin matrix for the L 2 -product using the basis functions 
(9i)i of X, namely (L4]x)f? = Jt A@i • Oj- Consequently, the notation [T + ] x + stands for 
the Galerkin matrix of the single layer operator defined on the triangulation X£. 



The EFIE can be discretized as follows : 



Find V + = (v!,.. .,v N +) T such that [T + ] x +V^ 

h 



-Id+K^x 



X7 

h 



and we have v + = J2iLi Viipf ■ The vector i? + v + is finally given by R + v + = YliLi Vifi- 
From now on, we denote by ( / the numerical computation of A^ described 



above. Using the former finite elements, the equation (15) takes the form of the linear 
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system 

where U = {ui,...,u n ) T and U = (u ,i, • • • , u 0tn ) T , with u = Y!,i=i u i ( Pi and u o = 

To precondition the DDM, we have mathematically proposed a multiplication by the 
operator T^. Numerically speaking, one wants to obtain a linear system close to the 
identity matrix. If we only multiplied the numerical vector by the matrix [T^jx^, we 
would obtain a linear system close to the mass matrix [Id]x h - Therefore, we have to do 
a numerical multiplication by the matrix [Id]^ [Ts]x h , in order to solve a system close 
to the identity matrix, and then better conditioned. As illustrated below, precondition- 
ing the method with the Galerkin matrix [Ts] x h is not enough to ensure an optimized 
convergence. One also needs to inverse the system by the mass matrix on the interface, 
which is realized through an iterative solution, of small numerical cost thanks to the 
sparsity of the matrix [Id]x h - This operation converts a vector whose coefficients are the 
L 2 -scalar products with the basis functions, into an amplitude vector (a vector whose 
coefficients are the coordinates in the basis functions) . 

We denote by DDM YO the original unpreconditioned equation related to the linear 
system (A J + A^ ) , 

DDM Yl is the equation with the left preconditioner being the Galerkin matrix of the 
single layer operator, 

[T s ]x h (( A J+) + ( A x-)) U = l T z]x h U . 

DDM Y2 is the equation with the left preconditioner being the Galerkin matrix of the 
single layer operator, with an additional inversion by the mass matrix [Id]x h , 

[Id]^[T E ] X/t + (A"-)) U = Mx h F*]x H U . 

DDM Y3 is the equation with the right preconditioner being the Galerkin matrix of the 
single layer operator, and an inversion by the mass matrix [Id]x h 

(( A I h +) + ( A x-)) Mxlmx h U = U . 

Let us remark that we have to solve two kinds of linear systems. The first one is the 
linear system arising from the DDM itself. The second one is made of the systems which 
come from the discretization of the integral equations inside each subdomain. In order 
to solve both of them, we use the GMRES algorithm. Notice that the numerical scheme 
is a doubly nested iterative method. 



14 



F. ALOUGES, J. BOURGUIGNON-MIREBEAU, AND D. P. LEVADOUX 



7.1. Validation of the method. First of all, we consider the degenerate case where 
there is no scattering object. We denote by E the sphere centered at the origin and of 
diameter lm, and we decompose M 3 into two subdomains: the interior and the exterior 
of S. Notice that Lemma |6.1| applies to this situation. The sphere meshed with 168 
DoF is shown on Fig. [5] (left), while the convergence curves are presented on Fig. [5] 
(right) at the frequency 68 MHz, and for the four equations above. Our first observation 
is that the unpreconditioned DDM YO converges really slowly in comparison with the 
three preconditioned DDM Yl, DDM Y2 and DDM Y3. DDM Y2 converges faster than 
DDM Yl, which lacks the inversion by the mass matrix. The results obtained with the 
right preconditioner (DDM Y3) are comparable with those of the left preconditioner 
(DDM Y2): both converge in as few as 4 iterations. 



DDM with analytic preconditioners 



Spherical mesh with 168 unknowns, frequency = 68 MHz, residue = le-6 




number of iterations 



Figure 5. Mesh (left) and convergence curves (right) to reach a residue 
of order 10~ 6 , for the artificial sphere at 68 MHz meshed with 168 DoF. 

In the next experiment, there is a scattering object which contains a cavity. This is 
the original setting intended for our study. First, we present the case of an object whose 
shape is close to a parallelepipedic box which is open at one of its extremities (Fig. 
[6]), and therefore exhibits a cavity. The interface S of this parallelepipedic box is a flat 
rectangle and is meshed with 102 DoF. For a residue of order 10 -6 , DDM Yl, DDM Y2 
and DDM Y3 converge respectively in 19, 13 and as few as 11 iterations, whereas the 
unpreconditioned method has not reached convergence after 1000 iterations. 
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DDM with analytic preconditioners 



Parallelepipedic mesh, frequency = 100 MHz, residue = le-6 




number of iterations 



Figure 6. Mesh (left) and convergence curves (right) to reach a residue 
of order 10~ 6 , for the parallelepipedic box at 100 MHz meshed with 102 
DoF on the interface X. 

On Fig. we compare the radar cross section (RCS) obtained by the four methods 
to the one obtained with the integral equation EFIE on the global mesh, without any 
artificial interface. Although this case is quite simple, there is no artefact due to the 
method. 

Curves of RCS for the mesh smallBox 

Frequency = 3200 MHz 



DDM Y0 

DDM Yl 

DDM Y2 




50 100 150 200 



angle of reception 



Figure 7. Curves of RCS for the mesh smallBox at a frequency of 3200 MHz. 
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7.2. Reliability of the method with respect to the frequency. Our third experi- 
ment illustrates the influence of the frequency increase on the convergence rate, for the 
sphere again, but with a finer mesh of the spherical interface E, of 3072 DoF. We choose 
to compare only DDM Y0 (without preconditioner) with DDM Y2 (with the left pre- 
conditioner). The number of iterations to reach a residue of order 10 -5 increases with 
the frequency for DDM Y0, whereas it remains stable (always 4 iterations) for DDM Y2 
(Fig. [8] and Tab. [TJ. Consequently, the convergence rate is not altered by the increase 
of the frequency, as illustrated on Fig. [5j on the right, and on Fig. [8j 



DDM with analytic preconditioner 

Spherical mesh with 3072 unknowns, at different frequencies, residue : 
„, 250 r 



le-5 




c 200 



150- 



100 



50 



G-© DDM Y0, unpreconditioned ' 
q-H DDM Y2, left preconditioner using [ld]~{-l}[T] 




H— B- 



100 



200 
frequency 



300 



400 



Figure 8. Mesh (left) and influence of the frequency increase (right) on 
the number of iterations to reach a residue of order 10~ 5 , for the artificial 
sphere meshed with 3072 DoF. 



Frequency (MHz) 


50 


68 


100 


150 


200 


250 


300 


360 


DDM Y0 


96 


101 


104 


160 


181 


189 


199 


216 


DDM Y2 


4 


4 


4 


4 


4 


4 


4 


4 



Table 1. Iterations count to reach a residue of order 10 5 depending on 
the frequency, for the sphere with 3072 DoF. 
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7.3. Reliability of the method with respect to the number of unknowns. 

7.3.1. Artificial spheres (no real object). 

We now refine the mesh of the sphere, passing from 3072 DoF to 5292 DoF. On this 
spherical mesh, at a frequency of 400 MHz, and to reach a residue of 10 -4 , DDM Y2 
converges in 28 iterations. The condition number of the linear system has obviously 
increased in comparison with the one of the former mesh (with 3072 DoF) (see Tab. [2]), 
leading to a smaller convergence rate. 

DDM with analytic preconditioner 



Spherical mesh with 5292 unknowns, frequency = 400 MHz, residue = le-4 




number of iterations 



Figure 9. Convergence curves to reach a residue of order 10 -4 , for the 
spherical mesh with 5292 DoF, at a frequency of 400 MHz, respectively 
for DDM Y0 (unpreconditioned) and for DDM Y2 (with analytic precon- 
ditioner). 



Equation 


Frequency 


Residue 


Number of unknowns 


Number of iterations 


DDM Y2 


360 MHz 




3072 


4 


DDM Y2 


400 MHz 


10" 4 


5292 


28 


Table 2. For the algorithm 


3DM Y2, comparison between the number 



of iterations needed to reach a given residue at a given frequency, respec- 
tively for the spherical mesh with 3072 DoF and for the spherical mesh 
with 5292 DoF. 

Nevertheless, looking at the convergence curves of the residues with respect to the 
iterations (Fig. [9]), we observe that the unpreconditioned DDM Y0 converges much 
slower than DDM Y2. In particular, after 28 iterations, DDM Y0 has not reached a 
residue of 5.10~ 2 , whereas DDM Y2 has reached a residue of 10~ 4 . As a conclusion, this 
preconditioner remains very efficient for finer geometries. 



18 



F. ALOUGES, J. BOURGUIGNON-MIREBEAU, AND D. P. LEVADOUX 



7.3.2. Hollow spheres (real objects). 



In this section, we illustrate the behavior of the algorithms (DDM YO to DDM Y3) 
when we refine the mesh of the scattering object. In that purpose, we consider a hollow 
sphere which constitutes the real scattering obstacle. The sphere is of radius 1 meter 
and is open for latitudes higher than 45 degrees, and is discretized with six different 



meshes of increasing precision (see Fig. 10 for the most refined mesh). 



The artificial interface needed for the DDM algorithm is chosen to be the missing cap 
of the sphere. Therefore, the interior and exterior problems consist in solving Maxwell 
equations inside and outside the sphere, respectively. They exchange data on the cap 
while, on the rest of the sphere, we have a Dirichlet type boundary condition. We give in 
Tab. [3] the number of unknowns respectively on the interface and on the whole sphere, 
for each considered mesh. Due to the size of the meshes, and contrarily to what has 
been done so far, we use a fast multipole method (FMM) to compress all involved linear 
systems. 



Name of the mesh 


Number of DoF 
on the interface (cap) 


Number of DoF on the spherical mesh 
of both subdomains 


hollowl2 


888 


5184 


hollow 15 


1380 


8100 


hollow20 


2440 


14400 


hollow25 


3800 


22500 


hollow30 


5460 


32400 


hollow35 


7420 


44100 



Table 3. Meshes of the considered hollow spheres. 
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Figure 10. Mesh (hollow35) of the hollow sphere of radius lm. The 
interface (not represented) possesses 7420 DoF while the whole sphere 
has 44100 DoF. 



Mesh : hollow35, interface : 7420 DoF, sphere : 44100 DoF 



Frequency = 400 MHz, residue = le-4 




20 40 

number of iterations 



Figure 11. Convergence curves to reach a residue of order 10 4 , for 
the mesh hollow35 of the hollow sphere, at a frequency of 400 MHz, 
respectively for DDM Y0 (unpreconditioned) and DDM Yl, DDM Y2, 
DDM Y3 (with analytic preconditioners). 
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Number of iterations with respect to the space discretization 

<.' Frequency = 400 MHz, residue = le-4 



3 40 

4— 

o 

01 

■v 

L_ 

(0 



30 



.20 



<u 10 



H— h DDM Yl, left preconditioner using [T] 

<^> DDM Y2, left preconditioner using [ld]~{-l}[T] 

A-A DDM Y3, right preconditioner using [ld]~{-l}=[T] 




1000 



2000 



3000 4000 5000 6000 
number of DoF on the interface 



7000 



8000 



Figure 12. Comparison of the number of iterations for the precon- 
ditioned equations DDM Yl, DDM Y2, DDM Y3, to reach a residue of 
10 -4 , at a constant frequency of 400 MHz. For each curve, each point cor- 
responds to one of the different meshes that we have considered, namely 
with 888, 1380, 2440, 3800, 5460 and 7420 DoF on the interface. 

We show in Fig. [TTJthe convergence rates for the four methods DDM Y0 to DDM Y3 
in the finest case (hollow35), for a frequency of 400 MHz. Once again, the unprecondi- 
tioned DDM Y0 converges much slower than the three other preconditioned equations 
(DDM Yl, DDM Y2, DDM Y3). For instance, DDM Yl (resp. DDM Y2, DDM Y3) 
reaches a residue of 10 -4 in 28 iterations (resp. 19, 24 iterations), whereas DDM Y0 has 
not yet reached a residue of 10~ 2 in 60 iterations. The explicit numbers of iterations for 
all meshes are provided in Tab. |4j 



Equation 
Mesh 


DDM Y0 


DDM Yl 


DDM Y2 


DDM Y3 


hollowl2 


152 


27 


17 


18 


hollowl5 


186 


26 


16 


16 


hollow20 


> 60 


27 


16 


19 


hollow25 


> 60 


27 


16 


19 


hollow30 


> 60 


28 


17 


22 


hollow35 


> 60 


28 


19 


24 



Table 4. Number of iterations to reach a residue of order 10 4 , at a 
constant frequency of 400 MHz, for the six meshes of the hollow sphere, 
for DDM Y0 (unpreconditioned) and DDM Yl, DDM Y2, DDM Y3 (with 
analytic preconditioners) . 
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Increasing the frequency of the problem to 1 GHz, only for the finest mesh (hollow35), 
does not deteriorate the method that much. Indeed, we show in Fig. [T3| and Tab. [5] 
that the three preconditioned methods DDM Yl, DDM Y2, DDM Y3 remain very 
competitive in comparison with DDM YO. 

Mesh : hollow35, interface : 7420 DoF, sphere : 44100 DoF 



Frequency = 1 GHz, residue = le-4 




c i i i i i i i i i d 

20 40 60 80 100 

number of iterations 



Figure 13. Convergence curves to reach a residue of order 10 -4 , for the 
mesh hollow35 of the hollow sphere, whose interface is meshed with 7420 
DoF, at a frequency of 1 GHz, respectively for DDM Y0 (unprecondi- 
tioned) and DDM Yl, DDM Y2, DDM Y3 (with analytic precondition- 
ers). 



Equation 
Mesh 


DDM Y0 


DDM Yl 


DDM Y2 


DDM Y3 


hollow35 


> 100 


59 


54 


53 



Table 5. Number of iterations to reach a residue of order 10 4 , at a 
frequency of 1 GHz, for the mesh hollow35, whose interface is meshed 
with 7420 DoF, for DDM Y0 (unpreconditioned) and DDM Yl, DDM Y2, 
DDM Y3 (with analytic preconditioners) . 
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8. Conclusion 

We have proposed a domain decomposition method associated with an efficient pre- 
conditioner, based on the restriction of the single layer operator on the interface between 
the subdomains. In each subdomain, the EFIE is solved at each iteration. The numer- 
ical results illustrate the very good behavior of the resulting preconditioned algorithm, 
which converges much faster than without preconditioning. 

Nevertheless, although the proposed method seems very encouraging, several difficul- 
ties need still to be overcome in order to make the method usable in real applications. 
First, the present formulation is restricted to the case where the EFIE is solved in each 
subdomain. Clearly, there is an obvious obstruction for resonant frequencies. In order 
to circumvent this issue, we have to generalize the approach for other formulations (e.g. 
CFIE, or the recent very efficient GCSIE methods [T], [2]). 

Another improvement direction consists in changing the coupling condition on the 
surface £ between the subdomains. For instance, when one takes impedant coupling 
boundary conditions, it is well known that the underlying problems are well-posed for 
any frequency [8]. Again, the GCSIE formalism, originally developed for metallic prob- 
lems, has recently been extended to impedant ones in [23] and |27j . and could prove to 
be very efficient. 

We plan to investigate those issues and even combinations of them in the foreseeing 
future. 

Acknowledgements. We would like to address special thanks to Jean-Marie Mirebeau 
for helpful suggestions. 
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